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Abstract 

In our recent work [H. Zhang, F.X. Trias, A. Oliva, D. Yang, Y. Tan, Y. 
Sheng. PIBM: Particulate immersed boundary method for fluid-particle in¬ 
teraction problems. Powder Technology. 272(2015), 1-13.], a particulate im¬ 
mersed boundary method (PIBM) for simulating fluid-particle multiphase flow 
was proposed and assessed in both two- and three-dimensional applications. 
In this study, the PIBM was extended to solve thermal interaction problems 
between spherical particles and fluid. The Lattice Boltzmann Method (LBM) 
was adopted to solve the fluid flow and temperature fields, the PIBM was re¬ 
sponsible for the no-slip velocity and temperature boundary conditions at the 
particle surface, and the kinematics and trajectory of the solid particles were 
evaluated by the Discrete Element Method (DEM). Four case studies were im¬ 
plemented to demonstrate the capability of the current coupling scheme. Firstly, 
numerical simulation of natural convection in a two-dimensional square cavity 
with an isothermal concentric annulus was carried out for verification purpose. 
The current results were found to have good agreements with previous refer¬ 
ences. Then, sedimentation of two- and three-dimensional isothermal particles 
in fluid was numerically studied, respectively. The instantaneous temperature 
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distribution in the cavity was captured. The effect of the thermal buoyancy on 
particle behaviors was discussed. Finally, sedimentation of three-dimensional 
thermosensitive particles in fluid was numerically investigated. Our results re¬ 
vealed that the LBM-PIBM-DEM is a promising scheme for the solution of 
complex fluid-particle interaction problems with heat transfer. 

Keywords: Lattice Boltzmann Method, Particulate Immersed Boundary 
Method, Discrete Element Method, Fluid-particle interaction, Heat transfer 


1. Introduction 


In recent years, numerical simulations based on combined Lattice Boltzmann 

a 


Method (LBM) pj and Discrete Element Method (DEM) [2j have gained pop¬ 
ularity in the fluid-particle interaction problems. The coupling scheme is quite 
attractive because the calculation of both sides is fairly local at the particle or 
even sub-particle scale. This feature provides more freedom on the particle ge¬ 
ometry and the choice of interaction laws. The commonest way for the coupling 
is to solve the fluid field at the particle scale while the solid particles are treated 
as moving boundaries [ 3 ]. The macro fluid quantities such as velocity and tem¬ 
perature are enforced to be coincident with the solid boundaries. This target 
is not easy to be achieved numerically since the solid boundary may not be 
at the same places with the lattice nodes. Therefore, the Immersed Boundary 
Method (IBM) Q is adopted to polish the stepwise representation and over¬ 
come the numerical oscillation when the moving boundary crosses the LBM 
grids. Meanwhile, the hydrodynamic force and torque exerted on the particle 
can be evaluated through the numerical correction. Finally, a proper particle 
tracking technique like the DEM Q is required to make the calculation cycle 
closed. 

For calculating the fluid-solid interaction force, there are several available 
schemes such as the penalty method [5], the direct forcing method [(| and the 
modified bounce-back rule proposed by Niu et al. Q. Comparing with the for¬ 
mer two, the third scheme is more efficient due to the fact that the forcing term 
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is simply evaluated based on the momentum exchange method. This scheme 


was thereafter tested in the Drafting-Kissing-Tumbling (DKT) problem [7] and 
simulation of several particulate systems |8[. However, neither of these two 
studies have considered more than two solid particles because the inter-particle 
collisions were roughly treated by the Lennard-Jones potential equation. In our 


previous work [9j, we reported a combined LBM-IBM-DEM method based on 


Q 


the scheme of Niu et al. [7|1. The new combined strategy was employed to simu¬ 
late the dynamic process of sedimentation of 504 particles in a two-dimensional 
square cavity. The advantage of this strategy is that the particle-particle in¬ 
teraction rules are governed by theoretical contact mechanics. Therefore, it 
has great potential to be a promising method since no artificial parameters 
are required during the calculation of both fluid-particle and particle-particle 
interaction forces. However, the major weakness of that is the low computa¬ 
tional efficiency. In the combined LBM-IBM-DEM method, one solid particle 
is represented by a set of small Lagrangian points. The interaction between the 
Lagrangian points and fluid lattice nodes provides detailed information of hy¬ 
drodynamic behaviors. On the other hand, the force and torque exerted on the 
particles is a summation action over all the relevant lattice nodes nearby. It is 


pointed out by Yu and Xu [lo| that the difficulty in particle-fluid flow modeling 
is mainly related to solid phase rather than fluid phase. A proper simplification 
on the fluid phase can be therefore tolerated especially when treating a system 
where inter-particle collisions dominate 


[ill . 12 . 13^ . For this reason, we further 


proposed a Particulate Immersed Boundary Method (PIBM) [14] to improve 
the computational efficiency of the original LBM-IBM-DEM scheme. The basic 
idea of the PIBM is to remove the constraints between the Lagrangian points 
and thus each Lagrangian point is treated as one single solid particle. Differ¬ 
ent to the conventional LBM-DEM based simulations [jjJ 16], the size of solid 
particles is allowed to be less than one fluid lattice in the PIBM. By doing so, 
detailed geometry of the solid particles is absent in the coupling whereas much 
more superior computational convenience is brought. The novel LBM-PIBM- 
DEM scheme has been successfully applied to three-dimensional simulation of 
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sedimentation process involving 8125 solid particles on a single CPU. In this 
study, we adopt the LBM-PIBM-DEM scheme to solve the thermal interactions 
between spherical particles and fluid. The LBM is adopted to solve the fluid 
flow and temperature fields, the PIBM is responsible for the no-slip velocity 
and temperature boundary conditions at the particle surface and the kinemat¬ 
ics and trajectory of the solid particles are evaluated by the DEM. To the best 
knowledge of the authors, no relevant research has been reported before. 

However, this is not the first attempt of using LBM or DEM to analyze 
the heat transfer phenomenon in a fluid-solid interaction system. Since He et 
al. [17] proposed a dual-LBM approach which uses a density distribution func¬ 


tion to simulate hydrodynamics meanwhile another temperature distribution 
function to simulate thermodynamics for heat transfer. The methodology has 
been widely adopted by various researchers (For example: Q0J [ 20 ! I 21 I and 


others). It is worthwhile mentioning that Han et al. [19[ proposed a numerical 
approach to account for the thermal contact resistance between contacting sur¬ 
faces of fluid and solid. They introduced a numerical case of a two-dimensional 
thermal cavity filled with solid particles to investigate the heat convection and 
conduction in the particulate system. The solid particles in the work of Han et 


al. Q 


were keeping stagnant and thus briefly played a role to construct com¬ 


plex solid boundary structures to test their model. Feng et al. 


22 


23] proposed 


a Discrete Thermal Element Method (DTEM) to model the heat transfer in the 
systems involving a large number of circular particles. Again, these work mainly 
focused on the heat conduction between solid-solid or solid-fluid whereas the dy¬ 
namic behavior of the solid particles were ignored. It is interesting to study the 
thermal convection of particulate flow in a fluid with intense inter-particle col- 

However, none of these work was conducted 
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based on LBM nor in three dimensional due to the enormous computational 
cost. In this study, we use the PIBM to conquer the limit. 

The remainder of the paper is organized as follows. To make this paper self- 
contained, the mathematics of the three-dimensional LBM, PIBM and DEM 
were briefly introduced in Sectional In Section [3j case studies of (1) Natural 
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Figure 1: Schematic diagram of the D3Q15 model [3 11 


convection in a two-dimensional square cavity with a concentric annulus, (2) 
Sedimentation of two-dimensional isothermal particles in fluid, (3) Sedimenta¬ 
tion of three-dimensional isothermal particles in fluid and (4) Sedimentation of 
three-dimensional thermosensitive particles in fluid were presented. The numer¬ 
ical results were discussed. Finally, conclusions were given in Section 0J 

2. Governing equations 

2.1. Lattice Boltzmann model with single-relaxation time collision 

We consider the simulation of the incompressible Newtonian fluids where 
the LBM-D2Q9 and D3Q15 models [l| are adopted for the two- and three- 
dimensional calculations, respectively. In this section, the equation systems of 
the D3Q15 model are presented. For the two-dimensional ones, the readers are 
referred to our previous study Q. The three-dimensional spatial distribution of 
the fluid velocities is shown in Figure [1] which can be expressed by 




(0,0,0)c 


a = 0 


e a =\ (±1,0,0)c, (0, ±1,0)c, (0,0, ±l)c a = 1-6 


( 1 ) 


(±l,±l,±l)c 


a = 7 - 14 


where c is termed by the lattice speed. The formulation of the lattice Bhatnagar- 




















Gross-Krook model is 


fa(r + e a 6 t ,t + 6 t ) = f a (r,t) - + F a S t (2) 

T f 

9air + e aStP + St) = g a (r,t) — ^ ^ - 9a ( i ) + Q a fi t (3) 

T g 

where f a (r, £) and g a (r, £) represent the fluid density and temperature distribu¬ 
tion functions, respectively, r = (x, y , z) stands for the space position vector, £ 
denotes time, and r/ and denote the non-dimensional relaxation times which 
can be calculated by 0 


\J Pvh c U c 

T/ “7Ia 


( 4 ) 


L c^C 


+ 0.5 


( 5 ) 


9 \fPrRac 2 5 t 

where c s = is the lattice speed of sound, L c and u c = \fgl3L c AT are the 
characteristic length and velocity, respectively, and Pr and Ra are the Prandtl 
and Rayleigh numbers, respectively. 


Pr = — (6) 

a f 

Ra = gpATL 
ajv f 

where Vf is the kinematic viscosity of fluid, af is the thermal diffusivity, g is the 
gravity, and /3 is the thermal expansion coefficient. The source terms, F a St and 
G a St , in Equations [2] and [3] are given in Section lT2l The equilibrium density 
and temperature distribution functions, /® g (r, t) and g^{r,t), can be written as 

fa (At) = PfU a [ 1 + 3(e„ • u) + ^(e a ■ u) 2 - ^ | u | 2 ] (8) 
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Figure 2: Schematic diagram of the PIBM. 


9 e a q (r,t) = Tu> a [l + 3(e a • u) + ^(e„ ■ u) 2 - ^ | u | 2 ] (9) 

where the value of weights are: uj o = 2/9, uj a = 1/9 for a = 1 — 6 and uj a = 1/72 

for a = 7 — 14. u denotes the macro velocity at each lattice node which can 

14 1 u c 

be calculated by u = ( / a e Q H—(AT( — ) 2 + Fs)St) / p/, the macro fluid 
a=o 2 c 

14 

density is pj = ^ f a and the macro temperature can be calculated by T = 

a=0 

14 1 

S 9a&a F ~QB^t' The discrepancy of the formulations of u and T here with the 
a=o 2 

conventional ones is due to the fact that they are modified by the momentum 

and heat flux Q. For example, Fb and Qb stand for the body force and heat 

u 

resp " tiveiy ’ ar< y 2 repre8e,ite ihe 

by temperature gradients [21]. 


2.2. PIBM 

It is worthwhile mentioning that the PIBMs for the no-slip velocity and 
temperature boundary conditions share a fairly similar logical manner. The 
common idea is to obtain an accurate expression of the velocity or temperature 
difference between the two phases and then use it to make further modification 
on the behavior of the fluid flow field and solid particles. In this study, the heat 
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conduction between the solid particles or particle and wall is not considered. 
For the sake of clarity, the two-dimensional schematic diagram of the PIBM is 
displayed with three-dimensional equation systems. As shown in Figure [2j the 
fluid is described using the Eulerian square lattices and the solid particles are 
denoted by the Lagrangian points moving in the flow field. The fluid density 
and temperature distribution functions on the solid particles are evaluated using 
the numerical extrapolation from the circumambient fluid points, 


f a (X l ,t) = L(X l ,r)-f a (r,t) (10) 

g«(X h t) = L{X u r) • g a (r,t) (11) 

where Xi(X,Y, Z) is the coordinates of the solid particles, the subscript T 
denotes the variable at the location of the Lagrangian particles. T(X/,r) is the 
three-dimensional polynomials, 

kmax 

Z Z{j n 

n=l,n\=k Zijk ~ Zijn 

( 12 ) 

where imax , jmax and k max are the maximum numbers of the Eulerian points 
used in the extrapolation as shown by blocks in Figure [2j With the movement 
of the solid particle, f a (Xi , t) will be further affected by the particle velocity, 

U p , 


L(X l ,r) = J2\ II 


X — Xijk 


n 


Vimk 


ijk 


l=l,l\=i Xi i k Xl i k ) \m=l,m\=j Vijk ViTnk 


fp(Xi,t + 6 t ) = f a (Xi,t) - 2uj a p f 


&aUp 


(13) 


and g a (. Xi , t) will be further affected by the temperature difference between the 


solid particle and fluid, AT 


I2l|, 


gp(Xi,t + 6 t ) = g a (Xi,t) - 2o; a AT^ (14) 

Ot 

in Equations [T3l and fT4l the subscript /3 represents the opposite direction of a, 
and h is the mesh spacing. Based on the momentum exchange between fluid and 
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particles, the force density, Ff(Xi,t), at each solid particle can be calculated 
using f a and fp, 


F f {X u t) = Y j ep[fp{X l ,t) - f a (X h t)\ (15) 

The effect on the flow fields from the solid boundary is the body force term F a S t 
in Equation [2j F a can be expressed by 

Fa = ^1 - u a ^3 —^2 - ~ ea ^j F z( r ’t) ( 16 ) 

where 

F B (r,t) = Ff (. Xi , t)Djjk (Tjjk — Xi)A p (17) 

i 

The heat source G a 5 t in Equation [3] is one dimensional and can thus be directly 
given as 


Go. — ^1 — (18) 

where 

Qb = ^==(^2AT— — Xi)A p ) (19) 

in Equations [l7| and [19j A p is the cross-sectional area of the particle which 
is given as A p = 0.257T d p , d p is the diameter of the particle. Dijj. c is used to 
restrict the feedback force to only take effect on the lattice nodes close to the 
solid particle and is given by 


DijkiTijk X{) — ^3 


%ijk Xl 

h 


S h 


Vijk 

h 


S h 


%ij k Z] 

h 


( 20 ) 


where 


Sh{a) 


\{l + co«s(^f)), when \ a |< 2 
0 , otherwise 


( 21 ) 
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Figure 3: The schematic diagram of (a) square cavity with a concentric annulus and (b) local 
zoomed view. 


On the other hand, the fluid-solid interaction force exerted on the solid particle 
can be obtained as the reaction force of Ff(Xi,t), 

Ffpi = —Ff(Xi,t)A p (22) 

2.3. Modeling of the particle-particle interactions 

Since heat transfer between particle-particle and particle-wall is not consid¬ 
ered, the remaining job is to monitor the trajectories of the solid particles and 
treat the inter-particle collisions properly. Based on the Newton’s second law 
of motion, the dynamic equations of the solid particle can be expressed as 


d 2 r 
771 dt? 

= (1- 

~~)d + Ff P i 

Pp 

(23) 

d 2 e 



(24) 

dt 2 

= T p 



where m and I are the mass and the moment of inertia of the particle, respec¬ 
tively. r is the particle position and 0 is the angular position, pf and p p are 
the densities of the fluid and particle, respectively. r p is the torque. Another 
considered force on the right hand side of Equation [23] is the fluid-particle in¬ 
teraction force Ff p i. When the particles collide directly with other particles or 
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Figure 4: Isothermal contours of temperature distribution in the square cavity at different Ra 
and Ar : (a)(top) Ra = 10 4 , Ar = 1.67, (b)(top) Ra = 10 4 , Ar = 2.5, (c)(top) Ra = 10 4 , 
Ar = 5; (b)(middle) Ra = 10 5 , Ar = 1.67, (b)(middle) Ra = 10 5 , Ar = 2.5, (b)(middle) 
Ra = 10 5 , Ar = 5; (c)(bottom) Ra = 10 6 , Ar = 1.67, (c)(bottom) Ra = 10 6 , Ar = 2.5, 
(c)(bottom) Ra = 10 6 , Ar = 5. 


the walls, DEM is employed to calculate the collision force. In this study, the 
particles and walls are directly specified by material properties in the simulation 
such as density, Young’s modulus and friction coefficient. When the collisions 
take place, the theory of Hertz 0 is used for modeling the force-displacement 
relationship while the theory of Mindlin and Deresiewicz 0 is employed for the 
tangential force-displacement calculations. For particle of radius Ri , Young’s 
modulus Ei and Poisson’s ratios z^, the normal force-displacement relationship 
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Figure 5: Isothermal contours of temperature distribution in the square cavity at different Ra 
and Ar : (a)(top) Ra = 10 4 , Ar = 1.67, (b)(top) Ra = 10 4 , Ar = 2.5, (c)(top) Ra = 10 4 , 
Ar = 5; (b)(middle) Ra = 10 5 , Ar = 1.67, (b)(middle) Ra = 10 5 , Ar = 2.5, (b)(middle) 
Ra = 10 5 , Ar = 5; (c)(bottom) Ra = 10 6 , Ar = 1.67, (c)(bottom) Ra = 10 6 , Ar = 2.5, 
(c)(bottom) Ra = 10 6 , Ar = 5. 


between the colliding particles reads 


F n = -E*R* 1/2 Sl /2 (25) 

3 

where the equivalent Young’s modulus and radius can be calculated by 1/E* = 
(1 — v\)/E\ + (1 — z %)/E 2 and 1/R* = 1/R\ + l/i? 2 , respectively. 

The incremental tangential force arising from an incremental tangential dis¬ 
placement depends on the loading history as well as the normal force 
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Ra 

Ar 

Present 

Hu et al. [20] 

Ren et al. [34] 

Moukalled and Acharya [35] 


5 

2.041 

2.038 

2.051 

2.071 

10 4 

2.5 

3.179 

3.184 

3.161 

3.331 


1.67 

5.213 

5.294 

5.303 

5.286 


5 

3.760 

3.778 

3.704 

3.825 

10 s 

2.5 

4.989 

4.917 

4.836 

5.080 


1.67 

6.193 

6.247 

6.171 

6.212 


5 

6.025 

6.095 

5.944 

6.107 

10 6 

2.5 

8.831 

8.934 

8.546 

9.374 


1.67 

11.857 

11.995 

11.857 

11.620 


Table 1: Comparison of computed average Nusselt numbers. 


AT = 8G*r a 9 k A5 t + (-l) k ^AF n (l - 0 k ) (26) 

where 1/G* = (1 — V\)/Gi + (1 — z/|)/G 2 , r a = \/S n R* is radius of the contact 
area. A St is the relative tangential incremental surface displacement, fi is the 
coefficient of friction, the value of k and Ok changes with the loading history. 


3. Results and discussions 


3.1. Natural convection in a two-dimensional square cavity with a concentric 
annulus 

Natural convection in a two-dimensional square cavity has been a popular 
benchmark case for the verification of one’s numerical tools on heat transfer 

, 2c|. In this subsection, the case is employed to 


through fluid materials 


35| 


m 


test the thermal LBM code coupling with the PIBM. An isothermal concentric 
annulus is planted in the center of the cavity as the heat source. The schematic 
diagram of the entire computational domain is shown in Figure [3^a). It should 
be stressed that the concentric annulus here is constructed by small isothermal 
particles with uniform size as shown in Figure [3fb). The positions of the solid 
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u o 1 
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Figure 6: Schematic diagram of the PIBM. 

particles are artificially set stagnant to prevent them following down under the 
action of gravity. The length and height of the cavity are 1, respectively. We 
consider three sizes of concentric annulus. The ratios of the cavity length to 
the annulus diameter Ar are 1.67, 2.5 and 5. We use 250 x 250 meshes to cover 
the whole computational domain. The diameter of the solid particles is h/2. 
Ra of interest is ranging from 10 4 to 10 6 and Pr = 0.71. The non-dimensional 
temperature is set 1 and 0 at the concentric annulus and the surrounding cold 
walls, respectively. The initial temperature of the stagnant fluid is 0.5. The 
remainder parameters relevant to the simulations are L c = 1 and u c = 0.25. 

The isothermal contours of temperature distribution in the cavity at different 
Ra and Ar are displayed in Figure 0J As shown, the temperature distribution 
is affected by both Ra and Ar. The strength of convection is stronger as Ra 
increases which gives rise to much more complex isothermal patterns. The ther¬ 
mal boundary at the lower surface of concentric annulus is thinner at the higher 
Ra for each Ar. As for the region just above the concentric annulus, a gradual 
changing of temperature can be observed at Ra = 10 4 which means that heat 
conduction mainly dominates at this region at low Ra. However, when Ra in¬ 
creases, a crown-like region with high temperature is formed due to the high 
effect of convection. Namely, the temperature distribution is more affected by 
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Figure 7: Positions of the 5000 particles without heat transfer at time (a) t = 2.5s, (b) 
t = 5.0s, (c) t = 7.5s, (d) t = 10.0s. 


the fluid velocity field at higher Ra. The streamlines in the cavity at different 
Ra and Ar are displayed in Figure [5j Similar to the temperature fields, the 
streamlines are symmetrically distributed about the perpendicular bisector of 
the concentric annulus. The influence of Ra on the flow field is obvious. As 
Ra increases, two more counter-rotating vortexes are shown between the upper 
surface of the concentric annulus and the top wall of the cavity when Ar = 1.67, 
the two vortexes besides the concentric annulus merge into a big one, respec¬ 
tively, when Ar = 2.5 and two more counter-rotating vortexes can be observed 
below the the concentric annulus when Ar = 5. All these observations have 
very good agreements with the previous numerical results under the same con- 
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Figure 8: Velocity distribution of the 5000 particles without heat transfer at time (a) t = 2.5s, 
(b) t = 5.0s, (c) t = 7.5s, (d) t = 10.0s. 


ditions 


m 


2D]. Quantitative comparison of computed average Nusselt numbers 


is shown in Table [TJ In this study, we calculate the average Nusselt number 
following Hu et al. 201 ]: 


_ i _ 

S t k(T hot - T^) 


(27) 


where k is the thermal conductivity. As can be seen from Table [TJ the LBM- 
PIBM scheme also does a good job of predicting reasonable Nu and thus is a 
promising method for simulating the free convection problems. 
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Figure 9: Velocity distribution of the 5000 particles with heat transfer at time (a) t = 2.5 s, 
(b) t = 5.0s, (c) t = 7.5s, (d) t = 10.0s. 


3.2. Sedimentation of two-dimensional isothermal particles in fluid 

From this subsection, fully coupled LBM-PIBM-DEM simulations are per¬ 
formed. The main target of this and next subsections are to investigate the 
temperature distribution feature in the cavity during particle sedimentation 
and the effect of the thermal buoyancy on particle behaviors. The current nu¬ 
merical results are compared with the ready-made cases in our previous study 
without considering heat transfer jl^ . 

In the two-dimensional simulation, we consider a 1 cm x 1 cm cavity with 
5000 two-dimensional particles. The calculating mesh for the LBM is 100 x 100. 
The diameter of the particles is 0.25 x 10 -2 cm or h/d p = 4. The initially spacial 
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(C) (d) 

Figure 10: Temperature distribution in the cavity at time (a) t = 2.5 s, (b) t = 5.0s, (c) 
t = 7.5s, (d) t = 10.0s. 


condition is exactly the same as the two-dimensional case in [l3| . Firstly, the 
5000 particles are randomly generated in the upper three-fifths domain and 
then deposit under the effect of the gravitational force, g = 9.8 m • s~ 2 . The 
density ratio between solid and fluid is 1.01, Ra = 10 4 , Pr = 0.71, L c — 1 and 
u c = 0.25. The non-dimensional temperature is set 1 and 0 at the solid particles 
and the four surrounding cold walls, respectively. The initial temperature of 
the stagnant fluid is 0.5. At last, the parameters responsible for the collision 
between the solid particles are E = 68.95 GPa and v = 0.33. For the sake 
of investigating the effect of the thermal buoyancy on particle behavior, two 
parallel simulations are carried out with and without considering heat transfer. 
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(c) 


(d) 


Figure 11: Velocity distribution of the 5000 particles with heat transfer at time (a) t = 12.5s, 
(b) t = 15.0s, (c) t = 17.5s, (d) t = 20.0s. 


In other words, the heat-introduced buoyancy would be ignored in the case 
without considering heat transfer. 

Figure [7| displays the particle distribution during the beginning 10.0s when 
heat transfer is not considered. As shown, the flow fluid at lower position is 
swallowed into the particle aggregate forming a fluid pocket with a mushroom 
shape. This typical two-dimensional phenomenon is regarded as the so-called 


Raylei 


^h-Taylor instability which has been well investigated by several stud- 
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13| . The solid particles prefer to settle around the fluid column 


until being shot up by the upthrust flow field. This process can be clearly seen 
from the distribution map of particle velocities as shown in Figure [51 Firstly, 
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Figure 12: Positions of the 8125 isothermal particles at time t = 0.0s. 

two vortexes with solid particles are formed in the interior of the cavity. Then, 
the height of these two vortexes rises meanwhile the downthrust of the particle 
streams introduces two more small vortexes close to the lower corners of the 
cavity due to the intimate interaction between the particle and fluid as shown 
in Figure |H(d). 

However, the sedimentation is fairly delayed when heat transfer is consid¬ 
ered. Figure [9] shows the velocity distribution of the 5000 particles with heat 
transfer at the same moments as Figure [8j Comparing with the case without 
heat transfer, the settling velocities of the thermal particles are significantly 
lower. This is due to the fact that the particle temperature is higher than the 
surrounding fluid. The fluid receives heat from the solid particles meanwhile 
the wall temperature is the lowest. Therefore, the fluid temperature increases 
at the location close to the solid particles whereas decreases close to the cold 
walls as shown in Figure [TUI It can be seen that the high temperature region is 
in conformity with the particle position. The temperature gradients make the 
density differences in the fluid occur. This gives rise to that the fluid in the cav¬ 
ity interior becomes less dense and rises. The solid particles at this region also 
rise with the ascending fluid which prevents the evolution in Figure 0 Besides 
the overall settling velocity, there is also large discrepancy on the particle dis- 
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(c) 


(d) 


Figure 13: Positions of the 8125 isothermal particles without heat transfer at time (a) t = 2.5s, 
(b) t = 5.0s, (c) t = 7.5s, (d) t = 10.0s. 


tribution patterns when heat transfer is considered. Instead of one fluid pocket 
with a mushroom shape in Figure 0 two branches are generated at the interface 
of solid and fluid as shown in Figure [9] The whole system is more unstable 
when heat transfer is introduced. 

Figure [Tl] displays the further evolution of the velocity distribution of the 
5000 particles when heat transfer is considered. As shown, the above-mentioned 
two branches in Figure [9]grow to two fluid pockets also with mushroom shapes. 
The two heads do not rise vertically any longer but towards the upper corners 
of the cavity. This phenomenon is not observed in the simulation without heat 


transfer before [36|, 151, M, LL3[. It seems that the thermal buoyancy totally breaks 
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(C) (d) 

Figure 14: Positions of the 8125 isothermal particles with heat transfer at time (a) t = 2.5 s, 
(b) t = 5.0s, (c) t = 7.5s, (d) t = 10.0s. 

the original rule and introduces more intensive fluid-particle interactions. 

3.3. Sedimentation of three-dimensional isothermal particles in fluid 

For the sake of conducting further investigation on the effect of thermal 
buoyancy on the particle behaviors, in this subsection, a three-dimensional 0.15 
cm x 0.15 cm x 0.15 cm cubic cavity is considered with six cold walls. The 
calculating mesh for the LBM isl5xl5xl5. The diameter of the solid particles 
is 0.5 x 10 -2 cm or h/d v = 2. The initial spacial set of the three-dimensional 

n 

case is the same as the three-dimensional case in [13] which also can be seen in 
Figure [12 8125 particles are regularly planted in the upper three-fifths domain 
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Figure 15: Minimum particle height versus time with and without heat transfer. 


of the cubic cavity with an identical separation distance, 0.001 cm. This leads 
to the local porosity, 0.719. The non-dimensional temperature is set 1 and 
0 at the solid particles and the six surrounding cold walls, respectively. The 
initial temperature of the stagnant fluid is 0.5. Other physical parameters of 
the fluid and particles are the same as the two-dimensional case in Section 13.21 
Similarly, two parallel simulations are carried out with and without considering 
heat transfer. 

Figure [13] displays the three-dimensional particle distribution during the 
beginning 10.0s when heat transfer is not considered. As shown, unlike the 
two-dimensional case, the three-dimensional particles in the center region settle 
more efficiently than others. The particles close to the walls move fairly slow as 
a result of the no-slip boundary condition employed between the fluid and wall. 
The particle matrix is rapidly hollowed out forming a downstream particle- 
constructed pestle. The particle-constructed pestle falls down directly until 
impacts on the bottom and then the particles scatter in all directions. This 
process has been reported in QH without considering heat transfer. 

Figure [IT] shows the particle distribution of the 8125 particles with heat 
transfer at the same moments as Figure [13] It can be seen that the particle 
behaviors are obviously different to the former case when heat transfer is not 
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(C) (d) 

Figure 16: Particle deposition velocity along the x-direction with and without heat transfer 
at time (a) t = 2.5 s, (b) t = 5.0s, (c) t = 7.5s, (d) t = 10.0s. 

considered. Instead of forming a particle-constructed pestle in the center region, 
the bottom surface of the isothermal particle aggregate is more fluctuant with 
bulges. This finding is in line with the two-dimensional one, the interface of 
the solid and lower fluid is more complex when heat transfer is considered. At 
t = 2.5s, the particle distributions exhibit completely contrary trends with and 
without heat transfer. It shows a hump in Figure [T3lai whereas a hollow in Fig¬ 
ure [2Ja). The discrepancy is due to the high thermal buoyancy in the interior, 
the particles in the corners have the priority to settle since the temperature in 
these regions is low. The subsequent sedimentation feature in Figure fl4l follows 
this trend where the bulging at the four corner region can be observed. 
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Figure 17: Isothermal surfaces in the cubic cavity at time (a) t = 2.5 s, (b) t = 5.0s, (c) 
t = 7.5s, (d) t = 10.0s. 


Figure [T5l shows the changing history of the minimum particle height in the 
two cases. This is a good measurement of the overall sedimentation efficiency. 
It can be seen that the two profiles are quite comparable before 2.4s. Then, 
the non-thermal particles begin to accelerate whereas the sedimentation of the 
thermal particles is significantly delayed by the buoyancy. Detailed comparison 
between the particle deposition velocity along the x-direction with and without 
heat transfer is given in Figure [T6l The contrary distribution of the particle 
deposition velocity validates the afore-mentioned observations. 

At last, the temperature distributions in the cubic cavity are displayed in 
Figure [l7| in terms of isothermal surfaces. The cavity is dissected at the bisector 


25 












(a) 




(c) (d) 

Figure 18: The 8125 thermosensitive particles colored by temperature at time (a) t = 0.0s, 
(b) t = 0.01s, (c) t = 0.03s, (d) t = 0.12s. 

of the x-direction which visualizes the temperature distribution feature inside. 
It is shown that the temperature at the interior of the particle aggregate is 
the highest. The thermal buoyancy at this region is also the highest and thus 
the solid particles possess the highest upwards velocities as shown in Figure HU 
The location of the heat core moves downwards with sedimentation of the solid 
particles. The temperature between the particle aggregate and surrounding 
walls changes gradually due to the low Ra adopted. 
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Figure 19: Temperature evolution histories of 7 particles. 


3.4 • Sedimentation of three-dimensional thermo sensitive particles in fluid 

In this subsection, three-dimensional thermosensitive particles are considered 
in the same model as Section ET31 Different to all the above cases where the 
temperature of the solid particles is constant. Here, the thermosensitive particles 
could lose or receive heat according to the surrounding temperature. The initial 
temperature of the 8125 thermosensitive particles is set randomly between 0 
and 1 as shown in Figure [18] (a). The walls are set cold with 0 temperature. 
The initial temperature of the stagnant fluid is 0.5. All the other computational 
parameters are the same as Section 13.31 

Since the walls are cold, it can be expected that the temperature of the 
whole system would reach a final steady-state which is the same as the wall 
temperature. This process can be clearly observed in Figure [18] which displays 
the particle temperature distribution at different time instants. Along with that 
the cold walls keep sucking heat, the temperature of the whole system drops 
rapidly to the steady-state value. Here, we monitor the temperature evolution 
histories of 7 solid particles randomly picked from the assembly. As shown in 
Figure [19j the different profiles again indicate the feature of the cooling process. 
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4. Concluding remarks 


In this paper, the LBM-PIBM-DEM scheme was employed to solve ther¬ 
mal interaction problems between spherical particles and fluid. The LBM was 
adopted to solve the fluid flow and temperature fields, the PIBM was responsi¬ 
ble for the no-slip velocity and temperature boundary conditions at the particle 
surface, and the kinematics and trajectory of the particles were evaluated by 
the DEM. 


Four case studies were implemented to certify the capability of the proposed 
coupling scheme. First, numerical simulation of natural convection in a two- 
dimensional square cavity with an isothermal concentric annulus was carried 
out for verification purpose. Then, sedimentation of two- and three-dimensional 
isothermal particles in fluid was numerically studied. The instantaneous tem¬ 
perature distribution in the cavity was captured. The effect of thermal buoy¬ 
ancy on the particle behaviors was discussed. Finally, sedimentation of three- 
dimensional thermosensitive particles in fluid was numerically investigated. Our 
results revealed that the thermal buoyancy has a great effect on the particle be¬ 
haviors both in two- and three-dimensional simulations. When heat transfer is 
considered, the interface between the solid particle aggregate and lower fluid is 
more unstable. The heat buoyancy could influence on the sedimentation effi¬ 
ciency of the solid particles very much. All the simulations demonstrate that 
the LBM-PIBM-DEM coupling scheme is a promising one for the solution of 
complex fluid-particle interaction problems with heat transfer. 

One critical assumption made in this study is that the heat conduction 
between solid particles or particle and wall is ignored. This assumption is 
reasonable when intense inter-particle collisions dominate the system, for ex¬ 
ample, the fore part of the sedimentation process where the collision process is 
instant. Nevertheless, this may be against those actual engineering pro cesses 
where the heat conduction between solids cannot be neglected 
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As 


shown in Section l3Tl the LBM-PIBM scheme is competent to pure natural con¬ 
vection problems. Further meshing on the solid particle can be performed to 
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calculate heat conduction through particles like Feng et al. 
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This work 


is interesting but will of course increase the calculation burden. An alternative 
approach is^to consider fluid flow and heat transfer just at the particle scale 
as done in 


26, 


27 


28 


Q. 


Therefore, it is in the authors’ opinion that the 


current LBM-PIBM-DEM scheme is more suitable for dynamic processes such 
as gas fluidization, hydrocyclone and pneumatic conveying. It can also be used 
to generate sub-particle information to support particle scale modeling. 
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